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A new method is presented for mesoscopic simulations of particle dispersions in nematic liquid 
crystal solvents. It allows efficient first-principle simulations of the dispersions involving many 
particles with many-body interactions mediated by the solvents. A simple demonstration is shown 
for the aggregation process of a two dimensional dispersion. 



PACS numbers: 61.30-Gz, 61.30. Jf, 61.20.Ja 



Dispersions of small particles in host fluids such as col- 
loidal suspensions and emulsions are of considerable tech- 
nological importance, and often appear in our everyday 
life in paints, foods, and drugs. Many kinds of exotic 
interactions between particles mediated by the host flu- 
ids are possible, including screened Coulombic 0, deple- 
tion [jl], fluctuation induced and surface induced || 
forces. A striking example occurs when spherical parti- 
cles are immersed in a liquid-crystal solvent in the ne- 
matic phase. For a single particle, the orientation of 
the solvent molecules is distorted due to the anchoring 
of the solvent molecules at the particle surface. Exten- 
sive studies have been done on this effect, and several 
characteristic configurations of the nematic field around 
a spherical particle have been identified @-|]. When the 
strength of anchoring is increased so that normal an- 
choring is preferred, the solvent changes from quadrupo- 
lar to dipolar symmetries around the particle. When 
more than two particles are immersed in the solvent, 
long-range anisotropic interactions are induced between 
particles due to elastic deformations of the nematic field 
fio| H. The anisotropic interactions can have a pro- 
nounced effect not only on the local correlations of the 
particles jl0|], but also on their phase behavior 14 IT]] 
and on their mechanical properties |l5|. 



Since analytical approaches for investigating these 
kinds of complex materials are extremely difficult, com- 
puter simulations are essential to investigate their static 
and dynamical properties. In most dispersions, the host 
fluid molecules are much smaller and move much faster 
than the dispersed particles. This enables us to assume 
that the host fluid is in local equilibrium for any given 
particle configurations, and thus it is usually a good idea 
to use some coarse grained mesoscopic descriptions for 
the host fluids rather than treating them as fully micro- 
scopic molecules [M. In the case of charged colloidal 
suspensions, a mesoscopic method for the first principle 
simulations can be derived by treating the counter ions 
as a charge density Jl9| . For the particle dispersions in 
nematic solvents considered here, the mesoscopic coarse 
grained free energy for the nematic solvent is well known 
to be the Frank free energy M, and the total free en- 



ergy T of the system consists of the following two parts; 
the bulk term T e j which presents elastic energy of the 
nematic and the surface term T s which determines an- 
choring of the nematic field at the particle surface. Let 
n(r) be the director, a common direction on which sol- 
vent molecules are aligned on average with a constraint 
n(r)| = 1. T can be given by functional of [n(r)] for a 
given particle configuration {Ri ■ • - Rat}; 



.f([n(r)];{R 1 -..R N }) 
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J dr [(V • n) 2 + (V x n) 2 ] + y j dS[l - (n • v)% (1) 



where K is the Frank constant with the single elastic 
constant approximation, W is the surface anchoring con- 
stant, and v is the unit vector normal to the colloidal 
surface HQ . The saddle-splay elastic term is not con- 
sidered in Eq. ([!]). The integral in the first term, T e i, 
runs over the whole solvent volume excluding the par- 
ticles, and that in the second term, J- s , runs over all 
solvent-particle interfaces. A simple scaling argument 
tells us T e i oc Ka d ~ 2 and T s oc Wa d ~ x with a and d 
being the particle radius and the system dimension, re- 
spectively, thus the physics should be determined by the 
ratio T sj J- e i cx Wa/K. Although this type of free en- 
ergy functional is sufficient for performing Monte Carlo 
simulations where only values of J- are needed for a given 
particle configurations, it is not useful for molecular dy- 
namics (MD) or Brownian type simulations because the 
coupling between solvent and the particles is given im- 
plicitly by limiting the integration space in both T e i and 
T s . This produces mathematical singularities at the in- 
terface when one calculates the force, ff s = —dF/dHi, 
acting on each particle mediated by the nematic solvents. 
Calculating the force is crucial for performing efficient 
simulations of many particle systems. Another serious 
problem of this type of functional is that in order to give 
correct boundary conditions at the particle-solvent in- 
terface, one has to use appropriate coordinates for per- 
forming grid-based numerical simulations rather than the 
usual Cartesian coordinates. This is generally difficult for 
particles with non-spherical shapes or for systems involv- 
ing many particles even when each particle has a sphcri- 
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cal shape. Also this makes the use of the periodic bound- 
ary condition difficult. 

To overcome these problems, we have modified Eq.(Q) 
by using a smooth interface between the solvent and the 
particles so that the coupling is given explicitly in the 
integrants through the interface. The new free energy 
functional we propose is 

T ([q(r)}-{R 1 ---R N }) 



ff s ({R 1 ---R N }) 
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where a,(3,j £ x,y,z and the summation convention is 
used. The explicit form of the interfacial profile fa be- 
tween dispersed particles and solvents is given by 



<M r ) = \ ftanh- ~ 
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under the director constraint (n a (r)) 2 = 1. One can 
perform this by numerical iterations such as the steep- 
est descent or the conjugate gradient method, ii) Once 
[QaP ( r )l ^ s obtained, the force acting on each particle me- 
diated by the nematic solvents follows directly from the 
Hellmann-Feynman theorem, 
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This form is very convenient because we can compute 
both dfa/dRi and 9(V ' a fa) / '9R at any time since fa is 
an analytical function of Rj. iii) Finally, we update the 

, .particle positions according to appropriate equations of 

^ 'motion such as 



with the particle radius a and the interface thickness 
f. Note that this reduces to Eq.(@) if Rc,£, -> 0. 
Very recently, a similar idea of using smooth interface 
was proposed for treating the hydrodynamic forces act- 
ing on particles dispersed in simple liquids [ pl| . In 
our case, the free energy is given by functionals of a 
traceless and symmetric second-rank tensor q a p(r) = 
n a (r)np(r) — 8 a p/d rather than the director n(r) to 
take into account the symmetry of the nematic director 
+n <-> — n automatically. The semi-empirical functional 
form l/R 2 tanh[i? 2 • • •] is applied in Eq.(|2|) to avoid the 
mathematical divergence of the elastic free energy density 
at the defect centers and to limit its value to A/ ~ K/ R 2 , 
which is the correct energy density difference between 
isotropic and nematic states, in the defect core regions of 
size R c . Another way to avoid the divergence would be to 
use the Landau-de Gennes type free energy with an or- 
der parameter Q a p(r) = Q(r)q a p(r), but this requires a 
prohibitively small lattice spacing near the defect points 

ii- 

The simulation procedure is as follows, i) For a given 
particle configuration {Ri • • • Rat}, we obtain the inter- 
face profile fa(r) by Eq.(||). Then we can calculate the 
stable (or meta-stable) nematic configurations [<Z^(r)] 
which satisfy the equilibrium condition 
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where ff p is the force due to direct particle-particle in- 
teractions (hard or soft sphere for instance), f$ and f/ 5 
are the hydrodynamic and random forces. Repeating the 
steps i)~iii) enables us to perform first-principles meso- 
scopic simulations for the dispersions containing many 
particles without neglecting many-body interactions. 

We have performed simple simulations for a two di- 
mensional (2D) system to demonstrate our simulation 
procedure. The system has 100 x 100 lattice sites in a 
square box with a linear length L = 100. Other physi- 
cal parameters are chosen rather arbitrarily as Rc = 1, 
a = 5, and £ = 2, where the unit of length is the lattice 
spacing I. Since the nematic configurations in 2D can be 
expressed by a single scalar field [0(r)], the tilt angle of 
the director against the horizontal (x-) direction, Eq.(|J) 
then reduces to 
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with q xx = cos 2 9 — 1/2, q yy = sin 9 — 1/2, and q xy = 
qyx = cos 9 sin 9. The boundary condition is fixed at 
9(r) — at the edge of the box to avoid rotations of the 
reference frame. We first calculated stable nematic con- 
figurations around a single particle for different Wa/K, 
and found two stable configurations. The first config- 
uration, which we refer to as weak anchoring, contains 
no topological defect. In the second configuration, which 
we refer to as strong anchoring, the particle is accompa- 
nied by two —1/2 charge point defects. Typical examples 
of the weak and strong anchoring are shown in Fig. 1(a) 
with Wa/K = 2 and Fig. 1(b) with Wa/K = 4, respec- 
tively. The distance between the defects and the particle 
center is about 1.3a. We note that both configurations 
possess quadrupolar symmetries, and the later would cor- 
respond to the Saturn ring configuration in three dimen- 
sional (3D) systems. Although in principle particles can 
be accompanied by one —1 charge hedgehog defect in 
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2D as well as in 3D, such configurations are unstable in 
the present 2D system since the elastic penalty of hav- 
ing m point defects with charge c scales as mKc 2 . This 
was directly confirmed by recent simulations with perfect 
normal anchoring [ p2| and also by our simulations. The 
total free energies T jW are plotted in Fig. 2 as functions 
of Wa/K. While both configurations can coexist in the 
narrow transition regime 2.2 < Wa/K < 2.5, our model 
predicts a clear first-order transition from the weak an- 
choring to the strong anchoring around Wa/K ~ 2.3. 

We next simulated the aggregation and ordering pro- 
cess of 30 colloidal particles after the isotropic to ne- 
matic transition of the solvent occurred. Here we used 
the periodic boundary condition and set Wa/K — 4 so 
that each particle is accompanied by two —1/2 charge 
defects. Other parameters are the same as in the previ- 
ous single particle case. The simulation was performed 
starting from a random particle configuration which is a 
typical configuration when the solvent is in the isotropic 
phase (K = 0). We then set K = 1 and calculated ff s 
according to the present procedure. The particle configu- 
rations were updated by numerically solving the steepest 
descent equation, 



fPS 



fPP 
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which is obtained by simply substituting cPR^/di 2 = 0, 
ff = 0, and ff = -CdRi/dt in Eq.(|). C = 1 is 
a friction constant and thus the off-diagonal compo- 
nents of the hydrodynamic interaction was not consid- 
ered. Here we obtain ff p = —dEpp/dHi from the 
repulsive part of the Lennard- Jones potential, Epp = 
"(2a/|r, - r,|) 12 - (2a/|r, - r,!) 6 + 1/4 
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to 



truncated at the minimum distance |r^ 
avoid the particles overlapping each other within the 
core radius ~ a. Snapshots from the present simulation 
are shown in Fig. 3(a) for an aggregation stage, and in 
Fig. 3(b) at a later time, where the particles are forming 
ordered clusters due to the anisotropic attractions be- 
tween them. Note is added that only up to two particle 
simulations have been done so far and simulations of 
more than three particles would be extremely difficult or 
almost impossible by other methods ever proposed. 

In summary, we have developed an extremely power- 
ful simulation method to investigate particle dispersions 
interacting via anisotropic solvents. We proposed a free 
energy functional which is suitable for MD type simula- 
tions. The following modifications have been made to the 
usual Frank free energy functional, i) The free energy is 
given by a functional of a tensor q rather than a vector 
n to take into account symmetry of the nematic director 
+n <-» — n. ii) The coupling between the nematic solvent 
and particles at the interfaces is introduced explicitly 
through a smooth interface so that we can analytically 
calculate the force acting on each particle mediated by 



the host by taking derivatives of the free energy according 
to the particle positions, iii) The value of the free energy 
density is limited semi-empirically to avoid a mathemat- 
ical divergence in the defect centers. We have performed 
demonstrations for a 2D dispersion and confirmed that 
the method works well even when the system contains 
point defects. Applications of this method to 3D systems 
should have no theoretical difficulties, but require heav- 
ier computation. This should allow the simulation of the 
chaining of the particles caused by the possible dipolar 
symmetry of the nematic configurations around a single 
particle. Although we have shown only simple demon- 
strations of the method by performing simulations of the 
2D system in this letter, simulations with physically more 
interesting situations such as systems with non-circular 
particles, asymmetric particle pairs with different parti- 
cle size, or particles with non-normal anchoring as well 
as more realistic simulations in 3D systems are now un- 
derway. 
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(a) (b) 




FIG. 1. Director configurations around a single parti- 
cle for (a) weak anchoring case without defect obtained at 
Wa/K — 2 and (b) strong anchoring case accompanied with 
two —1/2 charge point defects indicated with the crosses ob- 
tained at Wa/K = 4. The white disks indicate the particles 
with radius a = 5. Only 9% of the total system is shown for 
display purpose. 
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FIG. 2. The total free energies for the single particle cases 
as functions of the strength of the anchoring constant Wa/K. 
Our model predicts a first-order transition from weak to 
strong anchoring around Wa/K ~ 2.3. 
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FIG. 3. The aggregation and ordering process of colloidal 
particles when the solvent exhibit the isotropic [K — 0) to nc- 
matic (K = 1) phase transition at t = 0. Snapshots (a) in an 
aggregation stage (t = 10) and (b) at a later time (t = 100), 
where ordering of the particles are observed. Each particle 
is accompanied by two —1/2 charge point defects. Darkness 
presents the value of q£ x . Black and white correspond to 
q% x — and 0.25, respectively. Those correspond also to 
9 — 0.25-7T, 0.75-7T and 6 = 0, 0.5n, it as shown in the gradation 
map. 
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